      program plot2rot
      useutilities_module
      implicitnone
      integer,Parameter::irho0max=20
      integerIA,JMax,IMax
      Parameter(IA=150,JMax=70)
      Parameter(IMax=2*IA-1)
      integerNrDAT,NzDAT,NtDAT
      Parameter(NrDAT=IA,NzDAT=4*IA,NtDAT=JMax)
      integerNx,Ny,Narg,Iargc
      REALX(NrDAT),Y(NzDAT),Phi(NrDAT,NzDAT),W(1)
      character*80ResFile
      Nx=NrDAT
      Ny=NzDAT
      write(*,*)' N=',NrDAT,'   NzDAT=',NzDAT,'   NtDAT=',NtDAT
      write(*,*)' Nx=',Nx,'   Ny=',Ny
      callmypause
      write(*,*)' Arguments needed: phi*.res '
      Narg=Iargc()
      GOTO(09999,09998),Narg+1
      write(*,*)' Extra arguments.'
      Stop16
      GOTO09997
09999 CONTINUE
      ResFile='phirho0000.res'
      GOTO09997
09998 CONTINUE
      callGetArg(1,ResFile)
09997 CONTINUE
      write(*,'(3a)')'  ResFile=',ResFile(1:len_trim(ResFile)),' OK?'
      read(*,*)
      CALLFILLARY(ResFile,Nx,Ny,X,Y,Phi)
      CALLPLOT(X,Nx,Y,Ny,Phi,NrDAT,NzDAT,W,1.5,4,'Rcyl-axis','Z-axis','potential phi(Rcyl,Z)')
      end program plot2rot

      SUBROUTINE FILLARY(ResFile,Nx,Ny,X,Y,Phi)
      use utilities_module
      implicitnone
      integer,Parameter::irho0max=20
      integerIA,JMax,IMax
      Parameter(IA=150,JMax=70)
      Parameter(IMax=2*IA-1)
      integerNrDAT,NzDAT,NtDAT
      Parameter(NrDAT=IA,NzDAT=4*IA,NtDAT=JMax)
      character*80ResFile
      REALX(Nx),Y(Ny),Phi(Nx,Ny)
      integerNDIMP,NDP,NIP,Nx,Ny,ir,jt,i,j
      doubleprecisionpi,r(0:NrDAT),thet(0:NtDAT),Phirthet(0:NrDAT,0:NtDAT)
      doubleprecisionXd(Nx),Yd(Ny),radius,theta,PhiIntr0,PhiIntr1,Phid
      Parameter(pi=3.1415926535897932d+00)
      INTEGERM,LIWRK,LWRK
      PARAMETER(LIWRK=NrDAT+1+2*(NrDAT-2)*(NtDAT-2),LWRK=(NrDAT+7)*(NtDAT+7))
      INTEGERIFAIL,MX,MY,PX,PY
      DOUBLEPRECISIONC((NrDAT+1)*(NtDAT+1)),F((NrDAT+1)*(NtDAT+1)),FF(Ny+1),LAMDA(NrDAT+5),MU(NtDAT+5),WRK(LWRK),Xrad(0:Ny),Ythet(0:
     *Ny)
      INTEGERIWRK(LIWRK)
      EXTERNALE01DAF,E02DEF
      INTRINSICMAX,MIN
      write(*,'(3a)')' in _fillInput: ResFile=',ResFile(1:len_trim(ResFile)),' OK?'
      read(*,*)
      open(11,file=ResFile(1:len_trim(ResFile)))
      doir=0,NrDAT
      dojt=0,NtDAT
      read(11,*,end=10)r(ir),thet(jt),Phirthet(ir,jt)
      thet(jt)=min(thet(jt),pi/2.d0)
      enddo
      enddo
10    write(*,*)' after read ir=',ir
      write(*,*)' last ir, r(ir)=',NrDAT,r(NrDAT)
      write(*,*)' last jt, thet(jt)=',NtDAT,thet(NtDAT)
      callmypause
      WRITE(*,*)'E01DAF Program Results'
      MX=NrDAT
      MY=NtDAT
      dojt=0,MY
      doir=0,MX
      F((MY+1)*ir+jt+1)=Phirthet(ir,jt)
      enddo
      enddo
      IFAIL=0
      CALLE01DAF(MX+1,MY+1,r,thet,F,PX,PY,LAMDA,MU,C,WRK,IFAIL)
      WRITE(*,*)'               I   Knot LAMDA(I)      J     Knot MU(J)'
      doJ=4,MAX(PX,PY)-3
      IF(J.LE.PX-3.AND.J.LE.PY-3)THEN
      WRITE(*,99997)J,LAMDA(J),J,MU(J)
      ELSEIF(J.LE.PX-3)THEN
      WRITE(*,99997)J,LAMDA(J)
      ELSEIF(J.LE.PY-3)THEN
      WRITE(*,99996)J,MU(J)
      ENDIF
      enddo
      WRITE(*,*)'The B-Spline coefficients:'
      WRITE(*,99999)(C(I),I=1,MX*MY)
99999 FORMAT(1X,8F9.4)
99998 FORMAT(F5.2)
99997 FORMAT(1X,I16,F12.4,I11,F12.4)
99996 FORMAT(1X,I39,F12.4)
      doi=1,Nx
      Xd(i)=r(i)
      X(i)=Xd(i)
      enddo
      doj=1,Ny
      Yd(j)=(r(Nx)/real(Ny))*real(j)
      Y(j)=Yd(j)
      enddo
      doi=1,Nx
      doj=1,Ny
      radius=sqrt(Xd(i)**2+Yd(j)**2)
      theta=atan(Xd(i)/Yd(j))
      Xrad(j)=max(r(1),min(radius,r(NrDAT)))
      Ythet(j)=max(thet(1),min(theta,thet(NtDAT)))
      enddo
      M=Ny+1
      CALLE02DEF(M,PX,PY,Xrad,Ythet,LAMDA,MU,C,FF,WRK,IWRK,IFAIL)
      doj=1,Ny
      Phi(i,j)=FF(j)
      enddo
      enddo
      return
      end SUBROUTINE FILLARY
